Strong asymmetry for surface modes in nonlinear lattices with long-range coupling 
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We analyze the formation of localized surface modes on a nonlinear cubic waveguide array in the 
presence of exponentially-decreasing long-range interactions. We find that the long-range coupling 
induces a strong asymmetry between the focusing and defocusing cases for the topology of the surface 
modes and also for the minimum power needed to generate them. In particular, for the defocusing 
case, there is an upper power threshold for exciting staggered modes, which depends strongly on 
the long-range coupling strength. The power threshold for dynamical excitation of surface modes 
increase (decrease) with the strength of long-range coupling for the focusing (defocusing) cases. 
These effects seem to be generic for discrete lattices with long-range interactions. 
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Arrays of coupled optical waveguides and periodic pho- 
tonic lattices constitute a current area of intense research 
activity, due to the rich physical phenomena that arise 
when combining discreteness, periodicity, nonlinearity 
and surface effects [1 . Besides the interest stemming 
from the creation and controlling of the propagation of 
light beams for their potential use in multiport switching 
and routing of signals for envisioned all-optical devices, 
"discrete optics" has also recently become one of the fa- 
vorite tools for direct observation of phenomena associ- 
ated with discrete, periodic media, such as Bloch oscil- 
lations [2], Anderson localization [3 , discrete breathers 
and solitons [4], to name a few. 

A substantial amount of work has been devoted to 
the case of weakly-coupled nonlinear waveguide arrays, 
where the mode overlap between neighboring guides is 
small, and the nonlinearity is strictly local. Recent exper- 
imental and theoretical work in realistic systems such as 
dipole-dipole interactions in Bose-Einstein condensates 
(BEC) [5] and discrete light localization in nematic liq- 
uid crystals [6] , has stimulated research into the effects of 
nonlocal effects. In general, nonlocal nonlinearity tends 
to stabilize several types of solitons, such as dark solitons 
in 3D dipolar BEC [7] , chirp-imprinted spatial solitons in 
nematic liquid crystals [6 , optical vortex solitons [8] , ro- 
tating dipole solitons [9] and azimuthons [10 . The ef- 
fect of long-range dispersive interactions on the other 
hand, has received comparatively less attention. The ef- 
fect of power-law dispersion on anharmonic chains [TT] . 
as well as the inclusion of second-order coupling in optical 
waveguide arrays [12j [13] , suggests the onset of bistable 
effects. Although at first sight, the inclusion of long- 
range coupling would seem to lead to an increase of the 
power level needed to excite a localized mode [12], there 
are also some counterintuitive results for the case of a 
single nonlinear (cubic) defocusing impurity. There, a 
small addition of coupling to second nearest-neighbors, 
actually decreases the power threshold for the generation 
of a localized mode [14] . 

On the other hand, surface states have attracted con- 



siderable attention of the community during the last five 
years. Unlike the case of fundamental bulk modes, where 
there is no minimum power to excite them, for one- 
dimensional surface states there is a power threshold for 
their excitation. When only nearest-neighbors interac- 
tions are considered, this power is independent of the 
sign of the nonlinearity [15] . 
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Figure 1: iV-sites waveguide array with long-range coupling. 



In this work we examine the formation and stability 
of localized surface modes in a nonlinear optical waveg- 
uide array with realistic-looking long-range coupling (see 
FigfTl). We find a striking asymmetry between the behav- 
ior of the focusing and defocusing cases, as the coupling 
range is varied. Contrary to what occurs in a focusing 
case, for a defocusing nonlinearity an increase in cou- 
pling range actually reduces the amount of power needed 
to generate a surface localized stationary mode. This 
counterintuitive result also holds for the dynamical exci- 
tation of the surface mode from a narrow input beam. In 
addition, we found an upper threshold for the excitation 
of staggered states, effect that could be experimentally 
observed in current zig-zag arrays [13] . 

Let us consider a finite array of single- mode, nonlin- 
ear (Kerr) optical waveguides including higher-order cou- 
pling among sites. In the coupled-modes framework, the 
system is described by a discrete non-linear Schrodinger 



(DNLS) equation: 



.du n 
dz 



^ V n ,mUm + l\Un\ 2 Un = 0, (1) 



where u n is the amplitude of the waveguide mode in the 
n-th waveguide, z is the propagation distance along the 
array, 7 is the nonlinear parameter, and the coefficient 
V n ^m is the coupling between the n-th and ra-th guides. 
To be consistent with coupled-mode approach, we will 
model V nirn as V n , m = Ve -0 ^-™!- 1 ), where V is the 
usual coupling coefficient to nearest-neighbors and a > 
is the strength for the long-range interaction. A large a- 
value implies interaction with essentially one site (DNLS 
limit), while an small a increases the coupling range. The 
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Figure 2: (a) P versus A diagrams (including linear bands) 
for a = 6 (black), 2 (gray) and 1 (light-gray), for focusing 
and defocusing cases, (b) and (c) Linear modes for k — n and 
k = 0, respectively, (d) and (e) Nonlinear surface modes for 
points marked in (a). V = |^y| = 1 and N — 100. Black (gray) 
points in profiles denote u n > (< 0). 

power, defined as P = J^ n |n n | 2 , is a conserved quantity 
of model (fTl and we will use it to characterize nonlin- 
ear modes. We look for stationary solutions of the form 
u n( z ) — u n ex.p(iXz) of model (fiT), obtaining: 



Xu n = J2 Fe- a(|n - m| - 



'7<4 , 



(2) 



where u n G 3? and A is the propagation constant. To 
obtain the dispersion relation for linear plane waves, we 
set 7 = and insert a solution u n = Usin(kn) in pi), 
getting 



A(fc,a) 



/ e a cos k 



1 



cosh a — cos k 



(3) 



where k is the transversal wave number. Fig. J2|a) shows 
linear band regions for different values of a. The edges 
of these bands are located at Xmin = A(7r, a) = —2V/[l-\- 



exp(-a)] and X max = A(0,a) = 2V/[1 - exp(-a)], 
where the limit TV — >> 00 has been assumed. The width 
A maa; — A m i n = 4V/[1 — exp(— 2a)] increases as soon 
as coupling beyond nearest-neighbors is considered. As 
a consequence, the existence region for staggered solu- 
tions [A G {— 00, \mi n }] increases with a while the cor- 



responding region for unstaggered ones [A G {A^ 



>}] 



decreases. Figures [2^b) and (c) show profiles for A„ 
and X max with a staggered [( — l) n u n ] and an unstaggered 
(u n > Vn) topologies, respectively. 

Next, we compute nonlinear stationary surface solu- 
tions for focusing (7 > 0) and defocusing (7 < 0) non- 
linearities by implementing a multi-dimensional Newton- 
Raphson method [15 . A linear stability analysis reveals 
that the Vakhitov-Kolokolov criterion still holds in the 
presence of long-range coupling, i.e., dP/dX > implies 
stability. The P vs A curves for these modes show an 
important asymmetry between the focusing and defocus- 
ing cases [see Figj2ja)]: In the short-range coupling case 
(a = 6, black curves), power thresholds (Pth) for posi- 
tive and negative 7 are equal, like in a DNLS lattice [15]. 
However, when long-range coupling is relevant (a C 6), 
Pth increases as a decreases, in the focusing case. On 
the contrary, for 7 < this threshold decreases when a 
decreases [see gray and light-gray curves in FigJ2|a)]. 

An explanation for the P t h asymmetry can be the fol- 
lowing: We start from a surface profile, like the ones 
sketched in FigJ2|d) and (e), with the general form: 
tx {l, €,&£,...} with 1 > |e| > |0| > |£| > .... By 
inserting this ansatz in ([2| for n = 1, we get: A = 
V(e + fle~ a + £e~ 2a + ...) + juq. Since discrete soli- 
tons exist outside of linear bands, fundamental localized 
solutions would - in principle - bifurcate exactly at the 
frontiers of these bands, depending on the sign of non- 
linearity. Let us discuss first the unstaggered case and 
try to get an estimate for P t h in terms of a. It is well 
known that when the solution approaches the linear band 
(A — >• X max ), its power decreases and it becomes more and 
more extended (delocalized) [UH]. This implies that (in 
such a limit) e, /3, £, ... — >• 1 (this limit is exactly the op- 
posite to the one occurring for high level of power, where 
solutions are extremely localized and e, /?, £, ... —> 0). 



1/(1- 



lm- 



Therefore, 1 + e~ a + e~ 2a + e~ 3a + , 
plying that A — »> V/(l — e~ a ) -\-juq. However, this would 
imply that for uo — »■ 0, A < X maxi which is certainly a 
contradiction because the fundamental unstaggered so- 
lution could originate from the top of the band but not 
inside of it. As a consequence, at least 7^ ~ V/(l—e~ a ). 
Since power is directly proportional to Uq, we obtain the 
estimate P t h ~ 1/(1 — e~ a ). Thus, for 7 > 0, Pth 
will be a decreasing function of a, diverging at a = 
and remaining finite at a ^> 0. On the other hand, 
for 7 < 0, the situation is quite different. First of all, 
there is no a trivial transformation between unstaggered 
and staggered solutions as in the nearest-neighbor DNLS 
model. However, again, while the localized solution ap- 



proaches the linear band (A —> A min ), its power decreases 
and it becomes more delocalized, but now the solution 
is staggered. That implies a sign difference between 
nearest-neighbor amplitudes in the same way as the fun- 
damental linear mode located at Xmin [see FigMb)]. 
Therefore e, — /?, £, ... — »■ — 1. Now, we solve the sum: 
1 - e~ a + e~ 2a - e~ 3a + ... = 1/(1 + e~ a ), implying that 
A -> -V/(l + e~ a ) - \j\u%. For *i -► 0, A > A min , i.e 
a contradiction. Again, at least \j\uq « V/(l + e _Q; ), so 
Pt/i ~ 1/(1 + e _Q; ). Thus, for 7 < 0, P t h is an increasing 
function of a with a minimum at a = 0. Our analyti- 
cal estimates agree perfectly with the numerical behavior 
presented in Figj2|a). 




Figure 3: P versus A diagrams for modes centered at n = 1 
[dashed line, (a)], n — 2 [thick line,(b)], and n — n c [thin 
line,(c)]. V = \*y\ = a = 1 and N = 2n c = 100. 

We also computed localized solutions centered below 
the surface in order to detect the onset of the bulk phe- 
nomenology [see Fig. [3]. For 7 > 0, the power as a 
function of A shows the onset of a bistable curve for 
a < 1.69. This feature was observed before in the con- 
text of a zig-zag model [12j [13] , and seems to reflect an 
increase in effective dimensionality as soon as coupling 
beyond nearest-neighbors becomes important. The most 
salient feature is that in this case the threshold power 
to create a mode, behaves in a manner opposite to the 
usual DNLS. For example, for a = 1, Fig [3] shows that 
for 7 > 0, the minimum required power (Pth) f° r creating 
an unstaggered localized solution increases as the mode 
center is located away from the surface. In that sense, 
the system favors the localization of energy at the bound- 
ary for 7 > 0, contrary to the usual ID DNLS model [T5] 
(around a « 1.3, the DNLS phenomenology transforms 
into the long-range one). On the other hand, the system 
asymmetry is manifest for 7 < 0; the P t h for exciting a 
staggered localized mode decreases as the mode center is 
pushed away from the surface. Now, the system does not 
favor the generation of discrete surface solitons, as in the 
usual DNLS. 

Fundamental nonlinear modes are unstaggered for 7 > 
and staggered for 7 < 0. As the power content of 



the mode is increased, we find that unstaggered modes 
retain their character, as expected for high power solu- 
tions. However, contrary to what is expected, for stag- 
gered modes, a new power threshold "P up " appears where 
the staggered topology of the mode is lost. Figure [4] 
shows the existence region for unstaggered-staggered so- 
lutions in a-X space. We see that unstaggered solutions 
exist from some minimum A value (surface threshold) up 
to infinite. On the contrary, for 7 < the staggered 
mode is well-defined from the surface threshold value - 
which depends weakly on a - up to a second A value 
(with power P up ), that increases monotonically with a. 
Close to, but, beyond this second threshold, the mode 
is no longer staggered because although it retains some 
oscillations of the mode phases, it does not preserve a full 
staggered topology. As A decreases, the alternating phase 
topology is lost altogether. The white thick line separat- 
ing the light-gray and black regions is a result of asking 
the solution if us < 0, as an indicator of the change of 
topology in the central region. FigEFa) shows a mode 
example where all lattice sites are negative excepting the 
first one, i.e the mode is - by definition - not staggered. 
In the numerical continuation there is no evidence of this 
change on topology [see "sta" in Figj2|a) for 7 < and 
a = 1]. Curves are monotonous and the only way to 
observe this phenomenology is by taking a close look of 
phase structure. 

We can use an strongly-localized mode approximation 
to give an explanation for this unexpected behavior oc- 
curring for 7 < 0. This approach is valid when the prop- 
agation constant is far from the linear band, where the 
mode can be approximated as {u n } = uq{1, e, /3,0, ...}, 
and 1 ^> |e| ^> \/3\. We concentrate the analysis in the 




Figure 4: Existence regions in a-X space. The light-gray area 
corresponds to unstaggered (A > 0) and staggered (A < 0) 
solutions. The gray area represents the linear band while the 
white area denotes no solutions. The black region denotes 
solutions with no well-defined topology. Insets show examples 
of mode profiles in several regions. 



parameter /3, as an indicator of the long-range interaction 
effect. If we insert this ansatz in Q and solve it for site 
n = 3, we obtain: /3 « V(e + e- a )/(\-julP 2 ). From the 
anticontinuous limit, we know that high-power solutions 
consist essentially of one excited amplitude plus some 
exponentially small tails. Therefore, as a first approxi- 
mation, |A| > \j\uq. If 7 > also A,e > [see FigJ2|a) 
and (e)] implying that f3 > for any a. This shows us 
that, for a focusing case, solutions preserve their phase 
in the whole range of parameters. On the contrary, when 
7 < also A, e < [see FigJ2|a) and (d)], therefore the 
sign of /3 will depend on the balance |e| — e~ a . For a fixed 
a, this balance will be always negative for high-power so- 
lutions because e — » (anticontinuous limit). Therefore, 
for a large a an upper power threshold is expected ap- 
pearing at high frequencies; for smaller a, this threshold 
is expected to occur closer to the band because, there, e 
is also larger. This agrees perfectly with the thick white 
line of Figji} 




Figure 5: Output power fraction in Po-a space. Figures (a) 
and (b) correspond to 7 < and 7 > 0. Dark and light regions 
denote / = and / = 1, respectively. Insets: dynamical 
propagation for a = 1 and Po = 4. 

We also looked into the effects of long-range cou- 
pling on the dynamical evolution of an initially local- 
ized input beam. We solved Eq.(l) numerically with 
initial condition u n (0) = V^^ n ,i, where Po is the in- 
put power. For a given a and Po, we computed the 
space-averaged fraction of power remaining at the initial 
waveguide /, after a longitudinal propagation distance: 
/ = (P Zmax)- 1 f* max lu^z^dz. Results for / as a 
function of a and Po, are shown in FigJ5j in the form of a 
density plot. Its most striking feature is the diametrically 
opposite selftrapping behavior between 7 > and 7 < 0. 
While in the first case [FigJ5jb)], an increase in coupling 
range (decrease a) increases the threshold power for self- 
trapping, in the second case, a greater coupling range 
implies a smaller power threshold. This counterintuitive 
asymmetry becomes particularly strong around a ~ 1 
(see insets). These results are in complete agreement 
with the ones obtained for the stationary modes. 



Finally, we repeated all of the above studies on a sim- 
pler, but related model that is amenable to direct experi- 
mental probing: the zig-zag model [12 , and have verified 
the strong asymmetry effects for the formation of local- 
ized surface modes at both, the stationary and dynamics 
level. This opens the door to a direct experimental veri- 
fication of these effects. 

In conclusion, we have examined the formation of local- 
ized surface modes on a nonlinear waveguide array in the 
presence of realistic long-range interactions, and found a 
strong asymmetry between the focusing and defocusing 
cases for the mode topology and the minimum power to 
effect a localized surface mode. We believe these effects 
are generic to discrete nonlinear systems with long-range 
coupling. 
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